%This file calculates the equilibrium under untargeted subsidies in the full
%dynamic model for mature markets.



sl=1;

%ver=1;  %number of repetition
rng(ver);
nslot=20; %number of slots
ubound=0.15;  %calculated as 70% of the mean of ALL sunk cost (from new structual estimates), then divided by 20
yuma=(sl-1)*(ubound/nslot);


%### UPDATED on 12/08/2020: import market variables
textFileName_mktshift=sprintf('mktshift.csv');
MKTshift=dlmread(textFileName_mktshift);  %matrix: hsa year hhi totpop65 below50cat below75cat below100cat
MKTshift=sortrows(MKTshift,[1,2]);

  
%theta=[0 	-2.5349822	-4.8269864	-4.4024069 ...
%    	-3.6922329	-3.9341136 ...
%    0  	0.001956796 	0.002580179	  0.002569186 ...
%    	0.002515425	 0.002523947 ...
%    0.014364786 	0.16294816	0.1281231	0.17950069	-2.7699276	-2.9129867	-2.3893224	-2.0053693];   %the last row of parameters: mktdom per quartile (4 3 2 1), swc per quartile (1 2 3 4)


%theta=[0 -5.3121358	-5.1785215	-5.005499	-4.2737086	-4.4147539	-5.89177  ... %sunk costs
%   -4.5701719	-4.9060883	-3.3257925	-3.4972736	-4.0880094	-2.104362 ...
%    0 0.000751773	0.001047346	-0.000355248	-0.000534966	0.001187616	0.001131827 ... %variable costs
%    0.001270276	0.000107759	0.001047491	0.001054745	0.000681774	0.000649214 ...
%    0 -0.062058365	0.053012938	0.1467144	0.15329105	0.017731213	0.099901979 ... %vhhi
%    -0.056695029	0.10367065	0.10208316	-0.058478683	0.10910567	-0.011117688 ...
%    0 -0.003623524	-0.003845967	0.012251225	0.005590158	-0.014513229	-0.003046456 ... %vtotpop65
%   -0.018091252	0.005027301	-0.012294086	-0.010969303	0.0020008	-0.016308426 ...
%    0.19225046	0.53175016	2.5727366 ... %mktFEs
%    0.11259918	0.082341973	0.095907245	0.059739912	-2.0056423	-2.2808073	-1.9523567	-1.5906295];   %the last row of parameters: mktdom per quartile (4321), swc per quartile (1234)


%theta=[0  -5.1785215	-5.005499	-4.2737086  ... %sunk costs (pick vendor [3, 4, 5, 10, 12] from choice set {1:13}
%  	-3.3257925	-4.0880094	 ...
%    0  	0.001047346	-0.000355248	-0.000534966	... %variable costs
%    	0.001047491		0.000681774	 ...
%    0	0.053012938	0.1467144	0.15329105	 ... %vhhi    
%    	0.10208316	0.10910567	 ...        
%    0	-0.003845967	0.012251225	0.005590158	 ... %vtotpop65
%   	-0.012294086	0.0020008	 ...
%    0.19225046	0.53175016	2.5727366 ... %mktFEs
%    0.11259918	0.082341973	0.095907245	0.059739912	-2.0056423	-2.2808073	-1.9523567	-1.5906295];   %the last row of parameters: mktdom per quartile (4321), swc per quartile (1234)


%theta=[0 -5.1468019	-5.1489597	-4.650292	-4.1155525	-4.7014061	-5.7323459  ... %sunk costs
%          -4.6147978	-4.8042801	-3.306517	-3.633083	-4.1170821	-2.1809814 ...
%    0   0.000719732	0.001043895	-0.000342455	-0.000668024	0.001164287	0.001086266 ... %variable costs
%        0.001287248	2.95E-05	0.000975556	0.00100976	0.000719007	0.000728778 ...
%    0   -0.001570194	-0.019543905	0.12554481	0.13000603	-0.02172851	0.034346125 ... %vhhi
%        -0.18950474	0.072945226	0.04098723	-0.099828571	0.083168622	-0.026999884 ...
%    0   -0.005668858	-0.003234705	0.010732908	0.00571687	-0.011971641	-0.001144294 ... %vtotpop65
%        -0.020471505	0.00541441	-0.009380065	-0.008352656	0.002491854	-0.016889892 ...
%        0.18490543	0.52920777	2.6005711 ... %mktFEs
%        0.09900857	0.081858337	0.10076739	0.066720511	-1.8314467	-2.1632882	-1.896554	-1.5483258];   %the last row of parameters: mktdom per quartile (4321), swc per quartile (1234)

    
theta=[0 	-5.1489597	-4.650292	-4.1155525  ... %sunk costs
         -3.306517		-4.1170821	 ...
    0  	0.001043895	-0.000342455	-0.000668024	 ... %variable costs
       	0.000975556		0.000719007	 ...
    0  	-0.019543905	0.12554481	0.13000603	 ... %vhhi
        	0.04098723  	0.083168622	 ...
    0  	-0.003234705	0.010732908	0.00571687	 ... %vtotpop65
        	-0.009380065		0.002491854	 ...
        0.18490543	0.52920777	2.6005711 ... %mktFEs
        0.09900857	0.081858337	0.10076739	0.066720511	-1.8314467	-2.1632882	-1.896554	-1.5483258];   %the last row of parameters: mktdom per quartile (4321), swc per quartile (1234)


vhhi=theta(13:18);
vpop65=theta(19:24);


H=dlmread('ctft3fm2006_13chs.raw');  %matrix: hsa ahaid sysid bdtot newprechs (use the same set of hospitals except those HSAs with affiliated hospitals.
dropHSA=unique(H(H(:,3)>0,1));
for i=1:numel(dropHSA)
H(H(:,1)==dropHSA(i),:)=[];
end
H=sortrows(H,[1,4]);
data=H(:,[1,4]);


beta=0.95;
s=1:6;
ns=numel(s);
nh=3;
listmkt=unique(H(:,1));
nmkt=numel(unique(H(:,1)));


%{
bed=log(randi([20,250],nh*nmkt,1));
%bed=[4.02535169073515;5.29330482472449;3.25809653802148;4.41884060779660;4.64439089914137;5.26785815906333;4.81218435537242;3.97029191355212;4.59511985013459;5.33753807970132;4.74493212836325;4.88280192258637;5.41610040220442;5.11799381241676;3.36729582998647;5.15905529921453;4.54329478227000;3.66356164612965];
%bed=repmat([60;120;240],nmkt,1);
data=[kron((1:nmkt)',ones(nh,1)),bed];
data=sortrows(data,[1,2]);
%}
hid=1:nh;
[mktst2{1:nh}]=ndgrid(s); % the following 3 lines create a matrix with #rows=3^k and #columns=k and k=nh, the number of hospitals
mktst1=cat(1+nh,mktst2{:});
mktst=reshape(mktst1,[],nh);
clear mktst1 mktst2


%%%%%%%%%%%%%### UPDATED on 12/10/2020: find out initial states that are
%%%%%%%%%%%%%mature markets for each real market.
rng(1);
alldiff=zeros(size(mktst,1),1);
for i=1:size(mktst,1)
    onemkt=mktst(i,:);
    onemkt(onemkt==1)=[];
    if isempty(onemkt)==0 && numel(onemkt)==3
       onemkt=sort(onemkt); 
       if max(histc(onemkt,unique(onemkt)))==1
          alldiff(i)=1; 
       end
    end
end
tot_mature=find(sum((mktst>1),2)>2&alldiff>0);
select=datasample(tot_mature,nmkt,1,'Replace',true);
%id=randi([1,numel(tot_mature)],nmkt,1,'Replace',true);
mature_ini=[listmkt,select];


%%%%%%%%%%%%%### UPDATED on 04/16/2017: construct unconditional subsidies;
%%%%%%%%%%%%%conditional subsidies are constructed after mktdom is defined.
uncdsub=zeros(size(mktst,1),nh*ns);
cdsub=zeros(size(mktst,1),nh*ns);  %!!!defined after mktdom is defined
for i=1:nh
    onecol=(mktst(:,i)>1)*yuma;
    uncdsub(:,((i-1)*ns+1):(i*ns))=onecol(:,ones(ns,1));
end
%%%%%%%%%%%%%


eqmprob=zeros(size(mktst,1)*nmkt,ns*nh);  %### UPDATED on 01/10/2017: for eqm PROB 
eqmpi=zeros(size(mktst,1)*nmkt,ns*nh);  %### UPDATED on 05/08/2017: for eqm PROFIT 
for x=1:nmkt    
bd=data(data(:,1)==listmkt(x),2);
bdtype=(bd<=25&bd>0)+(bd<=80&bd>25)*2+(bd<=193&bd>80)*3+(bd>193)*4;
%### UPDATED on 12/08/2020: calculate the market shifter: sum of vhhi,
%vtotpop65, and mktFEs (MKTshift: hsa year hhi totpop65 below50cat
%below75cat below100cat)
hhi=MKTshift(MKTshift(:,1)==listmkt(x),3);
totpop65=MKTshift(MKTshift(:,1)==listmkt(x),4);
mktcat=MKTshift(MKTshift(:,1)==listmkt(x),5)*theta(25)+MKTshift(MKTshift(:,1)==listmkt(x),6)*theta(26)+MKTshift(MKTshift(:,1)==listmkt(x),7)*theta(27);
%mktcat2=MKTshift(MKTshift(:,1)==hsa,5)*theta(53);
%mktcat3=MKTshift(MKTshift(:,1)==hsa,6)*theta(54);
%mktcat4=MKTshift(MKTshift(:,1)==hsa,7)*theta(55);
%hsavar=vhhi(ones(period,1),:).*hhi(:,ones(ns,1))+vpop65(ones(period,1),:).*totpop65(:,ones(ns,1))+[zeros(period,1),mktcat2(:,ones(ns-1,1))]+[zeros(period,1),mktcat3(:,ones(ns-1,1))]+[zeros(period,1),mktcat4(:,ones(ns-1,1))];
%%%%%%%%%%%%
w=1/ns;
h1chs=w*ones(size(mktst,1),ns);
h2chs=w*ones(size(mktst,1),ns);
h3chs=w*ones(size(mktst,1),ns);
%h4chs=w*ones(size(mktst,1),ns);
for r=1:size(h1chs,1)
     if mktst(r,1)>1 
         h1chs(r,2:end)=1/(size(h1chs,2)-1);
         h1chs(r,1)=0;
     end
     if mktst(r,2)>1 
         h2chs(r,2:end)=1/(size(h2chs,2)-1);
         h2chs(r,1)=0;
     end
     if mktst(r,3)>1 
         h3chs(r,2:end)=1/(size(h3chs,2)-1);
         h3chs(r,1)=0;
     end
    % if mktst(r,4)>1 
    %     h4chs(r,2:end)=1/(size(h4chs,2)-1);
    %     h4chs(r,1)=0;
    % end    
end
hchs=[h1chs,h2chs,h3chs];
%%%%%%%%%%
%### UDPATED on 04/13/2017: change the transition probability to create
%incorrect belief for the focal hospital.
%{
wtp_h1=zeros(size(mktst,1),ns);
wtp_h2=zeros(size(mktst,1),ns);
wtp_h3=zeros(size(mktst,1),ns);
for r=1:size(mktst,1)
   wtp_h1(r,mktst(r,1))=1;
   wtp_h2(r,mktst(r,2))=1;
   wtp_h3(r,mktst(r,3))=1;
end
wrongtp=[wtp_h1,wtp_h2,wtp_h3];
%}
%%%%%%%%%%
v1=-mean(theta(1:ns))*1.5*rand(size(mktst,1),1);
v2=-mean(theta(1:ns))*1.5*rand(size(mktst,1),1);
v3=-mean(theta(1:ns))*1.5*rand(size(mktst,1),1);
%v4=-mean(theta(1:ns))*1.5*rand(size(mktst,1),1);
v=[mktst v1 v2 v3];  %excluding the hospital id
err=1;
listerr=zeros(10000,1);
newv=zeros(size(mktst,1),nh);
newp=zeros(size(mktst,1),ns*nh);
pi=zeros(size(mktst,1),ns*nh);   %### UPDATED on 05/08/2017: record the profit of each action in each market state for each hospital.
l=0;
while err>1e-6
l=l+1;
oldp=hchs;
oldv=v(:,4:6);
for i=1:nh
    %%%%%%%%%%
    %The following is to compute the net upfront payment.
    up=zeros(size(mktst,1),ns);
    ownchs=mktst(:,i);
    sw=(ownchs(:,ones(ns,1))~=s(ones(size(ownchs,1),1),:)).*s(ones(size(ownchs,1),1),:);
    up((sw(:)>0))=theta(sw(sw(:)>0)); 
    %cost=theta(1:ns);
    %up=cost(ones(size(mktst,1),1),:);
    %%%%%%%%%%
    %The following is to compute the switching cost.
    swc=zeros(size(mktst,1),ns);
    sw=(ownchs(:,ones(ns,1))~=s(ones(size(ownchs,1),1),:)).*(ownchs(:,ones(ns,1))>1);
    %swc(sw>0)=theta(end); 
    swc(sw>0)=theta(end-4+bdtype(i)); 
    swc(:,1)=0; 
    %swc=swc.*bd(i);
    %%%%%%%%%%
    %The following is to compute the interaction between the switching cost and large-hospital dummy.
    %swcNbig=(bd(i)>meanbd)*theta(end-2)*(swc~=0); 
    %%%%%%%%%%
    %The following is to compute the future cost-saving per bed.
    coef_bd=theta((ns+1):(2*ns));
    fvbd=bd(i)*coef_bd(ones(size(mktst,1),1),:);
    %%%%%%%%%%
    %The following is to compute the current value from bed size.
    coef_curbd=theta(ownchs+ns)';
    vcurbd=coef_curbd(:,ones(ns,1))*bd(i);
    %%%%%%%%%%
    %The following is to compute whether the current chosen is
    %market-leading.
    isdom=zeros(size(mktst,1),ns);
    curdom=zeros(size(mktst,1),1);
    curmsh=zeros(size(mktst,1),1);
    for ii=1:size(mktst,1)
        tmpmkt=mktst(ii,:);
        tmpmkt(i)=[];
        tmpmkt(tmpmkt==1)=[];
        if isempty(tmpmkt)==0
            occur=hist(tmpmkt,s);
            occur=reshape(occur,1,[]);
            mktshare=occur/sum(occur);
            curmsh(ii)=mktshare(ownchs(ii));
            mktdom=s(occur==max(occur));
            %### UPDATED on 12/10/2020: indicate whether a particular
            %vendor is market leading.
            s_trs=reshape(s,[],1);
            mktdom=reshape(mktdom,1,[]);
            isdom(ii,:)=max((~(s_trs(:,ones(numel(mktdom),1))-mktdom(ones(ns,1),:))),[],2)';                        
            curdom(ii)=max(mktdom-ownchs(ii)==0);
        end        
    end
    %vcurdom=theta(end-1)*curdom(:,ones(ns,1))*bd(i);
    %vcurdom=theta(end-1)*curdom(:,ones(ns,1));
    %vcurdom=theta(end-3-bdtype(i))*curdom(:,ones(ns,1));
    vcurmsh=theta(end-3-bdtype(i))*curmsh(:,ones(ns,1)); %### UPDATED on 09/03/2018: calculate the value of the current market share
    %### UPDATED on 12/08/2020: calculate the market shifters, depending on
    %the CURRENT choice
    mktval=vhhi(ownchs)'*hhi(1)+vpop65(ownchs)'*totpop65(1)+mktcat(1)*(ownchs>1);           
    %%%%%%%%%% ###UPDATED on 04/16/2017: define conditional subsidies.
    onecol=curdom*yuma;
    cdsub(:,((i-1)*ns+1):(i*ns))=onecol(:,ones(ns,1));
    %%%%%%%%%%
    %The following is to derive whether each choice is mkt_dom and the edv under each
    %possible future state.
    rhid=hid;
    rhid(rhid(:)==i)=[];
    tempmst=mktst;
    tempmst(:,i)=[];  
    rst=unique(tempmst,'rows');
    ddom=zeros(size(rst,1),ns);
    fv=zeros(size(rst,1),ns);
    fmid=[i rhid]; 
    for j=1:s(end)  
        fmst=[j*ones(size(rst,1),1) rst];
        for k=1:size(fmst,1)
            %The following is to compute the edv under each possible future state.
            onemst=[fmid;fmst(k,:)];
            nnmst=sortrows(onemst',1)';
            idx= v(:,1)==nnmst(2,1)&v(:,2)==nnmst(2,2)&v(:,3)==nnmst(2,3);
            fv(k,j)=v(idx,i+nh);
            %The following is to compute the ddom under each possible
            %future state.
            config=fmst(k,:);
            config(:,1)=[];   %### UPDATED on 12/14/2016: exclude the choice of the focal hospital
            config(config(:)==1)=[];
            if isempty(config)==0
               occur=histc(config,s)';
               occur=reshape(occur,1,[]); 
               mktshare=occur/sum(occur);  %### UPDATED on 09/03/2018: construct market share
               %occur=occur+affocr;
               mktdom=s(occur==max(occur));  %a row vector
               ddom(k,j)=max(mktdom-j==0);          
            end
        end
    end
    %%%%%%%%
    %The following is to compute the transition probability.
    temphchs=hchs;  %64*(4+4+4)
    c1=ns*(i-1)+1;
    c2=ns*i;
    temphchs(:,c1:c2)=[];
    tp=zeros(size(mktst,1),size(rst,1));
    for t=1:size(rst,1)
        tp(:,t)=temphchs(:,rst(t,1)).*temphchs(:,rst(t,2)+ns);  
    end
    %csv=up+swc+beta*fvbd+beta*tp*ddom*(theta(end-1))+beta*tp*fv;   %### UPDATED on 12/13/2016: reverse back to one logsum    
    %csv=up+swc+beta*tp*ddom*(theta(end-1))+beta*tp*fv;   %### UPDATED on 12/13/2016: reverse back to one logsum    
    %csv=up+swc+beta*fvbd+beta*tp*ddom*(theta(end-1))+beta*tp*fv;   %### UPDATED on 12/13/2016: reverse back to one logsum  
    %csv=vcurmsh+vcurbd+up+swc+beta*tp*fv;
    %csv=vcurmsh+vcurbd+up+swc+[zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))]+beta*tp*fv;   
    %csv=vcurmsh+vcurbd+up+swc+[zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))]+cdsub(:,((i-1)*ns+1):(i*ns))+beta*tp*fv; 
    %csv=vcurmsh+vcurbd+up+swc+[zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))]+beta*tp*fv-swc.*isdom; 
    csv=vcurmsh+vcurbd+up+[zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))]+beta*tp*fv;         
    %csv=vcurdom+vcurbd+up+swc+beta*tp*fv;
    %### UPDATED on 12/10/2016: for logsum case, also needs to adjust
    %available actions that can be taken according to the inital state,
    %"mktst".    
    expcsv=exp(csv);
    expcsv(v(:,i)>1,1)=0;
    %expcsv(expcsv(:,1)>0,1)=1;
    ecsv=log(sum(expcsv,2));
    %ecsv(v(:,i)==1,1)=0;
    %%%The following is to update ev by computing the prob-weighted csv.
    %ecsv=sum(hchs(:,c1:c2).*csv,2);
    %%%The following is to update ev by picking the max of csv.
    %csv(ownchs(:,1)~=1,1)=-Inf;
    %ecsv=max(csv,[],2);
    newv(:,i)=ecsv;
    updfv=zeros(size(rst,1),ns);
    for m=1:s(end)  
        fmst=[m*ones(size(rst,1),1) rst];
        for g=1:size(fmst,1)
            %The following is to compute the edv under each possible future state.
            onemst=[fmid;fmst(g,:)];
            nnmst=sortrows(onemst',1)';
            idx= v(:,1)==nnmst(2,1)&v(:,2)==nnmst(2,2)&v(:,3)==nnmst(2,3);
            updfv(g,m)=ecsv(idx,1);
        end
    end
    %updcsv=vcurmsh+vcurbd+up+swc+beta*tp*updfv;
    %updcsv=vcurmsh+vcurbd+up+swc+ [zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))] +cdsub(:,((i-1)*ns+1):(i*ns))+beta*tp*updfv-swc.*isdom;  
    updcsv=vcurmsh+vcurbd+up+ [zeros(size(mktst,1),1),mktval(:,ones(ns-1,1))] +beta*tp*updfv;      
    %updcsv=vcurdom+vcurbd+up+swc+beta*tp*updfv;
    %updcsv=up+swc+beta*tp*ddom*(theta(end-1))+beta*tp*updfv;
    %updcsv=curdom(:,ones(ns,1))*theta(end-1)+up+swc+beta*tp*updfv;    
    %updcsv=up+swc+swcNbig+beta*fvbd+beta*tp*(ddom*(theta(end-1)+yuma)+updfv);
    %updcsv=up+swc+swcNbig+beta*fvbd+beta*tp*(ddom*(theta(end-1)+yuma)+theta(end-2)*(bd(i)>meanbd)*ddom+updfv);
    expupdcsv=exp(updcsv);
    expupdcsv(v(:,i)>1,1)=0;
    %expupdcsv(expupdcsv(:,1)>0,1)=1;
    sumexpuv=sum(expupdcsv,2);
    ppp=expupdcsv./sumexpuv(:,ones(ns,1));
    %}
    hchs(:,c1:c2)=ppp;
    newp(:,c1:c2)=hchs(:,c1:c2);
    pi(:,c1:c2)=updcsv; 
end
   d1=abs(newp-oldp);
   d2=abs(newv-oldv);
   err=sum(d1(:))+sum(d2(:));
   v(:,4:6)=newv;
   listerr(l,1)=err;     
end
eqmprob((size(mktst,1)*(x-1)+1):(size(mktst,1)*x),:)=newp;
eqmpi((size(mktst,1)*(x-1)+1):(size(mktst,1)*x),:)=pi;
%alleqmpr(x,:)=newp(:,1);
end



%forward simulate with estimated equilibrium probability
period=20;
alleqmpr=eqmprob;
%Indicator for each hospital
allmkt_fchs=zeros(nmkt*nh,period);
allmkt_fdom=zeros(nmkt*nh,period);
allmkt_furfdom=zeros(nmkt*nh,period);
allmkt_fpi=zeros(nmkt*nh,period);
%Average across all hospitals in each market
avg_fadp=zeros(nmkt,period);
avg_fdom=zeros(nmkt,period);
avg_furfdom=zeros(nmkt,period);
avg_fpi=zeros(nmkt,period);
for x=1:nmkt
    %compute future choices
    eqmpr=alleqmpr((size(mktst,1)*(x-1)+1):(size(mktst,1)*x),:);
    %w=reshape(eqmpr(1,:)',ns,[])';
    w=reshape(eqmpr(mature_ini(x,2),:)',ns,[])';  %### UPDATED on 12/10/2020: forward simulation start from the initial state in mature markets    
    cc=cumsum(w,2);
    cc(:,end)=1;
    chs=s(sum(bsxfun(@ge,rand(nh,1),cc),2)+1);
    chs=reshape(chs,[],1);
    allmkt_fchs(((x-1)*nh+1):(x*nh),1)=chs;
    all_fadp(((x-1)*nh+1):(x*nh),1)=(chs>1);
    avg_fadp(x,1)=mean(chs>1);
    ftmpmkt=chs;
    ftmpmkt(ftmpmkt==1)=[];
    %compute mktdom
    if isempty(ftmpmkt)==0
       occur=hist(ftmpmkt,s);
       occur=reshape(occur,1,[]);
       mktdom=s(occur==max(occur));
       allmkt_fdom(((x-1)*nh+1):(x*nh),1)=max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2); 
       avg_fdom(x,1)=mean(max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2));
    end
    %compute further definition of mktdom
    if isempty(ftmpmkt)==0
       occur=hist(ftmpmkt,s);
       occur=reshape(occur,1,[]);
       mktdom=s(occur==max(occur));
       if max(occur)<2
           mktdom=0;
       end
       allmkt_furfdom(((x-1)*nh+1):(x*nh),1)=max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2); 
       avg_furfdom(x,1)=mean(max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2));
    end
    %compute profits
    ftmppi=eqmpi((size(mktst,1)*(x-1)+1):(size(mktst,1)*x),:);
    onemkt_ftmppi=ftmppi(1,:);
    onemkt_ftmppi=reshape(onemkt_ftmppi',ns,[])';
    allmkt_fpi(((x-1)*nh+1):(x*nh),1)=onemkt_ftmppi(sub2ind(size(onemkt_ftmppi),(1:size(onemkt_ftmppi,1))',chs));
    avg_fpi(x,1)=mean(onemkt_ftmppi(sub2ind(size(onemkt_ftmppi),(1:size(onemkt_ftmppi,1))',chs)));
    %compute the same things from period 2 onwards
    for pd=2:period
        future_idx= v(:,1)==chs(1)&v(:,2)==chs(2)&v(:,3)==chs(3);
        w=reshape(eqmpr(future_idx,:)',ns,[])';
        cc=cumsum(w,2);
        cc(:,end)=1;
        chs=s(sum(bsxfun(@ge,rand(nh,1),cc),2)+1)';
        chs=reshape(chs,[],1);
        allmkt_fchs(((x-1)*nh+1):(x*nh),pd)=chs;
        avg_fadp(x,pd)=mean(chs>1);
        ftmpmkt=chs;
        ftmpmkt(ftmpmkt==1)=[];
        if isempty(ftmpmkt)==0
           occur=hist(ftmpmkt,s);
           occur=reshape(occur,1,[]);
           mktdom=s(occur==max(occur));
           allmkt_fdom(((x-1)*nh+1):(x*nh),pd)=max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2);   
           avg_fdom(x,pd)=mean(max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2));
        end        
        if isempty(ftmpmkt)==0
           occur=hist(ftmpmkt,s);
           occur=reshape(occur,1,[]);
           mktdom=s(occur==max(occur));
           if max(occur)<2
               mktdom=0;
           end
           allmkt_furfdom(((x-1)*nh+1):(x*nh),pd)=max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2); 
           avg_furfdom(x,pd)=mean(max(~(chs(:,ones(size(mktdom,2),1))-mktdom(ones(size(chs,1),1),:)),[],2));
        end
        onemkt_ftmppi=ftmppi(future_idx,:);
        onemkt_ftmppi=reshape(onemkt_ftmppi',ns,[])';
        allmkt_fpi(((x-1)*nh+1):(x*nh),pd)=onemkt_ftmppi(sub2ind(size(onemkt_ftmppi),(1:size(onemkt_ftmppi,1))',chs));
        avg_fpi(x,pd)=mean(onemkt_ftmppi(sub2ind(size(onemkt_ftmppi),(1:size(onemkt_ftmppi,1))',chs)));
    end    
end
tot_cdsub=(allmkt_fdom>0)*yuma;
tot_uncdsub=(allmkt_fchs>1)*yuma;


output=[mean(avg_fadp,1);mean(avg_fdom,1);mean(avg_furfdom,1);mean(avg_fpi,1);sum(tot_uncdsub,1)];


eqmname=sprintf('eqm3Dsub_mature_pre.csv');
resultname=sprintf('result3Dsub_mature_pre_ver%d.csv',ver);


dlmwrite(resultname,output,'delimiter', ',', 'precision', 8)
dlmwrite(eqmname,eqmprob,'delimiter', ',', 'precision', 8)




